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When non-adsorbing polymers are added to an isotropic suspension of rod-like colloids, the col- 
loids effectively attract each other via depletion forces. We performed Monte Carlo simulations to 
study the phase diagram of such rod-polymer mixture. The colloidal rods were modelled as hard 
spherocylinders; the polymers were described as spheres of the same diameter as the rods. The 
polymers may overlap with no energy cost, while overlap of polymers and rods is forbidden. 

Large amounts of depletant cause phase separation of the mixture. We estimated the phase 
boundaries of isotropic-isotropic coexistence both, in the bulk and in confinement. To determine 
the phase boundaries we applied the grand canonical ensemble using successive umbrella sampling 
[J. Chem. Phys. 120, 10925 (2004)], and we performed a finite-size scaling analysis to estimate 
the location of the critical point. The results are compared with predictions of the free volume 
theory developed by Lekkerkerker and Stroobants [Nuovo Cimento D 16, 949 (1994)]. We also 
give estimates for the interfacial tension between the coexisting isotropic phases and analyse its 
power-law behaviour on approach of the critical point. 

PACS numbers: 82.70.Dd, 64.70.Ja, 68.35.Rh, 05.10.Ln 



I. INTRODUCTION 

Non-adsorbing polymers are often added to colloidal 
suspensions in order to modify the effective interactions 
between the colloids. By this means phase transitions 
can be driven, e. g., in order to facilitate the isolation 
of the colloids^ -. In particular, mixtures of viruses and 
polymers are widely used in experiments on colloidal liq- 
uid crystals (see, for instance, the recent review by Dogic 
and Fraden 9 -). 

As certain viruses, such as Tobacco Mosaic Virus and 
fd- virus, can to a first approximation be regarded as long 
cylinders, they are good model systems for liquid crystals. 
Currently there is much interest in the non-equilibrium 
behaviour of these systems - in particular in the effect 
of shea r 11 ! 12 ! 15 ! 16 . In order to interpret the experimental 
results, detailed theoretical knowledge on the equilibrium 
phase diagram is needed. Therefore we present in this ar- 
ticle a computer simulation study of the phase behaviour 
of rods and spheres, both in the bulk and in confinement. 

The Asakura-Oosawa-Vrij mode l 19 ! 29 is a useful model 
for mixtures of polymers and spherical colloids. The 
polymers are assumed to be freely interpenetrable with 
respect to each other, while there is a hard-core repulsion 
between the colloids as well as between the colloids and 
the polymers. As the interaction energy is then either 
zero or infinite, the phase behaviour is purely of entropic 
origin. 

In this article we discuss a similar model for mixtures 
of rod-like colloids and polymers: a mixture of hard sphe- 
rocylinders with length L and diameter D and freely in- 
terpenetrable spheres with the same diameter. 

This system was studied with liquid-state integral 
equation theory 5 -^ and free volume theory^ as well as 
with computer simulations 4 - ! 17 i 20 . In recent years vari- 



ous regimes of size-ratios and concentrations have been 
discussed, e. g. the packing properties at very high 
concentrations or the behaviour of small rods, which 
act as a depletant on large spheres. Here we focus on 
dispersions of rods and spheres of similar diameters at 
low concentrations. The phase diagram of this model 
was studied within free volume theory by Lekkerkerker 
and Stroobants 14 . Details of this approach will be de- 
scribed in section [TTJ Li and Ma recently computed the 
effective interaction between two rods by Monte Carlo 
simulation^. Bolhuis and coworkers 4 - as well as Savenko 
and Dijkstra^ determined the bulk rod-sphere phase di- 
agram by simulation (via thermodynamic integration). 
Both have given results for rods of aspect ratio L/D = 5 
and spheres of several diameter values. In order to 
avoid simulating the spheres explicitly, Savenko and Di- 
jkstra used an exact expression for the effective Hamilto- 
nian, which was numerically evaluated during the Monte 
Carlo simulation for each rod configuration. Bolhuis and 
coworkers 4 - modelled the spheres explicitly in their Gibbs 
Ensemble Monte Carlo simulations to study fluid-fluid 
separation, while for the other parts of the phase dia- 
gram they used an effective expression for the rod-rod 
interaction. 

Here we present results for the full rod-sphere sys- 
tem, which were obtained in the grand canonical en- 
semble. The successive umbrella sampling method^ was 
employed to determine the grand potential hypersurface 
of the system. This allowed us to predict the phase 
boundaries of isotropic-isotropic coexistence with much 
higher accuracy than the studies mentioned above, which 
used thermodynamic integration. In particular, we deter- 
mined the critical point by a finite-size scaling analysis. 
We also studied the influence of confinement on the phase 
behaviour. This is relevant for experiments under shear. 
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In section |TT] we briefly review free volume theory for 
the rod-sphere mixture. In section IIIII the simulation 
method is introduced. In section IIVI we show results 
for phase diagrams and interfacial tensions and compare 
them to the theoretical predictions. Section [V] contains 
a summary and a discussion. 



II. FREE VOLUME THEORY 

We briefly review the theoretical approach to rod- 
sphere mixtures which was introduced by Lekkerkerker 
and Stroobants 14 in 1994. The starting point is the ther- 
modynamic potential in the semi-grand canonical ensem- 
ble, where the number of rods is fixed, while the num- 
ber of "penetrable hard" spheres is set by the chemi- 
cal potential of spheres in a hypothetical reservoir that 
is held in equilibrium with the system. The potential 
f2 semi (AT r , V, T, /i s ) of such a mixture can be written as: 
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where F(N r , V, T) is the free energy of N r rods in a vol- 
ume V at temperature T, and the second term accounts 
for the perturbation due to the addition of N s spheres 
from a reservoir, where the chemical potential is kept at 
pf . Osmotic equilibrium requires the chemical poten- 
tials of the system and the reservoir to be equal. The 
chemical potential of an ideal gas of spheres is connected 
to the density p a via 



exp 



(2) 



Thus, the number of spheres in the system depends on 
pf. The only influence the rods have on the spheres is 
the reduction of the free volume: 



(3) 



where Vf ree is the volume accessible for spheres under the 
assumption that the rod configurations are undistorted 
upon adding spheres and p R is the number density of 
spheres in the reservoir. The free volume fraction a is 
defined as: 



a 



V 
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and it can be calculated within scaled particle theory^. 
Hence, the expression for the semi-grand canonical po- 
tential is reduced to: 



(5) 



The chemical potential as well as the osmotic pressure 
IL. of the rods in the mixture can be obtained from the 



thermodynamic relationships: 
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where v r = ttD 3 (2 + 3L/D)/12 is the volume of a sphe- 
rocylinder of length L and diameter D. 

The virial expansion of the free energy of a system 
of hard spherocylinders can be calculated using scaled 
particle theory 7 : 



= j~~t ( const + Hv) + *[f] + n 2 y + % 2 
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where r\ r = p r v r is the volume fraction of rods. 
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and t is the overall length-to-diamcter ratio of the sphe- 
rocylinders 



L + D 
D 
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In eqs. |8] and [TOl c[f] and p[f) are functionals of the 
orientational distribution function /: 

a[f] = /"/(u)ln[47r/(u)]rfu (12) 

Pif] = ^J ysin[<f>(u,u')]/(u)/(u')dudu' (13) 

where u is the unit vector along a particle's axis. 

The normalised orientational distribution function of 
a single rod / is adapted to minimise the semi-grand 
canonical potential. In the nematic case it can be 
obtained either numerically 23 or from the Gaussian 
approximation 1 ^, but here we are interested in the 
isotropic case, i. e., all orientations u are equally probable 
and the distribution is uniform: 



(14) 



Therefore, in our case, the functionals are reduced to: 



<r\f\ = 
P[f] = 1 



(15) 
(16) 
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The expression for the free volume fraction calculated 
with scaled particle theory reads: 



1 
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where, in the isotropic case, for "penetrable hard" 
spheres with diameter D: 
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Thus, the chemical potential and the osmotic pressure 
of the rods in the mixture are functions of the rod vol- 
ume fraction and of the chemical potential of spheres in 
the reservoir. At coexistence, they are equal. From this 
condition we have obtained the phase diagrams, which 
are compared to our simulation results in section ITVl 



III. MODEL AND SIMULATION METHOD 

In our Monte Carlo simulations we model the colloids 
as hard spherocylinders of length L and diameter D. The 
polymers are approximated as spheres of the same diam- 
eter, which are hard with respect to the rods and freely 
interpenetrable among each other. 

We performed simulations in the grand canonical en- 
semble, where the volume V and the chemical potentials 
fi r of the spherocylinders and /x s of the spheres are fixed. 
The temperature T is formally also fixed, but irrelevant, 
since it only sets the energy scale. Simulations were per- 
formed in a box with edges L z > L y = L x > 3L and 
periodic boundary conditions. In a rectangular box con- 
figurations at coexistence form preferably such that the 
interfaces are parallel to the small faces of the box (see 
figure [T|) . This simplifies the analysis of the interfacial 
tension. The finite-size effects were examined in a cubic 
box. For the study of confinement effects on our model 
system we choose a geometry where the box dimensions 
were L x — L y > 3L and L z — 3L with periodic boundary 
conditions in the x- and y-directions and hard walls in the 
z-direction. To speed up the simulations we employed a 
cell system for efficient overlap detection of anisotropic 
particles^,. In this approach the simulation box is cut 
into cubes of side lengths > D. Whether a rod intersects 
a cube or not can be computed very fast. The volume 
which needs to be checked for overlaps then contains at 
the most a few particles. 



A. Cluster move 

In order to efficiently equilibrate canonical simulations 
of spherical colloids, Biben and co-workers introduced a 




liquid-like ■ \j gas-like ■ liquid- 

FIG. 1: Configuration snapshot for L/D = 5 at coexistence. 



cluster movei. Later this move was extended to simula- 
tions in the grand canonical ensemble 2 ^. Here we intro- 
duce a version that is adjusted to the case of spherocylin- 
ders. 

If we attempt to insert a spherocylinder into a simu- 
lation box full of spheres, it will certainly overlap with 
some of them and the move will be rejected. The clus- 
ter move combines the insertion of a spherocylinder with 
the removal of the spheres, which it otherwise would over- 
lap with. Thus it increases the probability of accepting 
an insertion move. At the same time, the removal of a 
spherocylinder is combined with the insertion of a cer- 
tain amount of spheres in the void left by the removed 
spherocylinder . 

Care has to be taken of the acceptance probabilities in 
order to ensure detailed balance. Here, we only give the 
most relevant equations. The reader who is interested in 
details is referred to the work by Vink and Horbach 26 , 
who have explained the process thoroughly for spheri- 
cal colloids. For spherocylinders we have to account for 
an additional degree of freedom - the orientation. The 
Metropolis probability of accepting a removal of a sphe- 
rocylinder in a cluster move is: 



acc(N r -> JV r - 1) 



N r m(z s V dz ) n ' 
AttV z r n s \ 
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(19) 

where AE is the potential-energy difference between the 
initial and the final configuration - which is infinite for 
an overlapping configuration and zero otherwise. The 
volume of the depletion zone around a spherocylinder is 
V dz = 7r(2D) 3 (3L/(2D)+2)/12. The factors 1/(4tt) ineq. 
IT91 and 47r in eq. [2T] are due to the orientational degree 
of freedom. The fugacities z r / s of spherocylinders and 
spheres respectively are related to the chemical potentials 
via z r j s = exp(/3(x r / s ), where /3 = l/(ksT) is the inverse 
temperature. n s is the number of spheres to be inserted 
into the void. n s is drawn uniformly from the interval 
[0, to), where to is an integer given by to = 1 + max[l, 
int(z s Vd z + a^/ z s Vdz)}, with a a positive constant of order 
unity. 

If only one spherocylinder is removed and no spheres 
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are inserted, the acceptance probability is reduced to: 



acc{N r — * N r — 1) = min 
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The Metropolis probability of accepting an insertion of 
a spherocylinder in a cluster move is: 



acc(N r 



N r 



1) = min 
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if n s < m and otherwise. Here n s is the number of 
spheres which need to be removed to produce a void, in 
which the spherocylinder can be inserted. If only one 
spherocylinder is inserted and no spheres are removed, 
the acceptance probability is reduced to: 
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FIG. 2: Grand potential ln[P(r; r )] for spherocylinders of as- 
pect ratio L/D — 5 at sphere fugacities z 3 = 1.0125 (crosses), 
1.0625 (stars) and 1.175 (squares). Afl indicates the height 
of the grand potential barrier between the coexisting phases. 



B. Phase boundaries from grand-canonical 
simulations 

We determined the phase boundaries from the prob- 
ability distribution P(N r ), which is the probability to 
observe N r rods in the mixture for given chemical po- 
tentials jj, r and /i s . For a given value of fi s , we searched 
for the value of fi r at which the distribution is bimodal. 
However, the simulations need not be performed right 
at the coexistence chemical potential (which is usually 
not known beforehand). In fact, the chemical potential 
can be set to any value /J^" 1 and then renormalised to 
coexistence via 

ln[P^ (N r )] = ln[P^ to (JV r )] + ( Mr - /4 im )AV (23) 

such that the areas under the two peaks in P(N r ) are 
equate. We used the successive umbrella sampling 
method^ to determine P(N r ). This technique allowed 
us to sample regions between the two bulk phases where 
P(N r ) is very low. 

Figure [2] shows the logarithm of the probability dis- 
tribution P(j] r ) as a function of the rod volume fraction 
7] r = N r v r /V with v r = ttD 3 (2 + 3£/Z))/12, which - up 
to an additive constant - is equal to the grand potential 
of the system. The locations of the maxima of the peaks 
are the bulk volume fractions of rods at coexistence. For 
a large fugacity z s (squares), and hence a large concen- 
tration of spheres, there are two clearly separated peaks 
indicating a phase transition which is strongly of hrst 
order. With decreasing concentration of spheres (stars 
and crosses) the effective attraction between the rods be- 
comes weaker. Hence the peaks broaden and move closer 
until eventually the critical point is reached. 

Figure [1] shows a typical configuration at coexistence. 
Due to the periodic boundaries two slabs separated by 
two interfaces have formed. 

In figure one can also clearly see that an advanced 
biasing scheme such as successive umbrella sampling is 



necessary to bridge the huge differences in probabilities 
between the pure bulk states and the states in the two- 
phase coexistence region (which show up via a horizontal 
part of the f2 vs. r\ r curve, since a change of r\ r just 
amounts to a change of the sizes of the gas-like and liquid- 
like domains, without changing the areas of the interfaces 
between them). 



IV. RESULTS 
A. Bulk 

1. Phase diagram 

The phase diagram is presented choosing the fugacity 
of the spheres and the rod volume fraction as variables. 
The fugacity z s is related to the sphere reservoir volume 
fraction via rif = z s v s , where v s — irD 3 /6 is the volume 
of a sphere. Explicit implementation of the spheres allows 
an easy transformation into the frame of (r) s ,T] r ), where 
t] s is the sphere volume fraction in the system. In free 
volume theory the actual concentration of spheres follows 
as 7y s = a7]f , with a from eq. [T7J 

Figures [3] and H] show phase diagrams for mixtures of 
rods with aspect ratios L/D — 3 and 5 and spheres. The 
solid lines are free volume theory predictions. The sin- 
gle phase, the isotropic mixture of rods and spheres, is 
marked "I" . The region of the phase diagram, where the 
gas-like and liquid-like isotropic phases coexist is named 
"I+I" . The two almost vertical lines are theoretical pre- 
dictions for the phase boundaries of the isotropic-nematic 
transition^. 

Since free volume theory is based upon a mean-field 
approximation, and fluctuations, which are especially rel- 
evant near the critical point, are ignored, we expect its 
predictions to deviate from the simulation results there. 
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FIG. 3: Bulk phase diagram for a mixture of spherocylinders 
with aspect ratio L/D = 3 and spheres of diameter D. The 
solid lines are predictions of the free volume theory. The filled 
square marks the critical point. The inset shows the phase 
diagram in the (z s , ?) r )-plane. The almost vertical lines in the 
inset indicate the coexistence region of the isotropic-nematic 
transition. 
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FIG. 5: Comparison of the bulk phase diagram for a mixture 
of spherocylinders with aspect ratio L/D = 5 and spheres 
of diameter D to previous work. The open triangles are the 
results of the present work, the filled square marks the critical 
point. The open squares are the data obtained by Bolhuis and 
coworkers^. The solid lines are predictions of the free volume 
theory. 
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FIG. 4: Bulk phase diagram for a mixture of spherocylinders 
with aspect ratio L/D = 5 and spheres of diameter D. The 
solid lines are predictions of the free volume theory. The filled 
square marks the critical point. The inset shows the phase 
diagram in the (z a , ?) r )-plane.The almost vertical lines in the 
inset indicate the coexistence region of the isotropic-nematic 
transition. 



Away from the critical point the predictions of the free 
volume theory approach the simulation results. As anti- 
cipated, the theory underestimates the volume fraction of 
spheres in the liquid phase considerably. This is due to 
the fact that the depletion forces change the local struc- 
ture of the fluid2i - an effect which is not included in 
the calculations of the free volume accessible to spheres. 
On the gas-branch of the coexistence region, where the 
amount of rods is negligible, the theoretical predictions 
agree well with the simulation results. 

Figure [5] shows a comparison of our results with the 
data obtained in previous computer simulations on the 



fluid-fluid separation by Bolhuis and coworkers^ (open 
squares). The errorbars of their data are Aij r ~ 0.02. 
Thus, the results do not coincide within the errorbars. 
We attribute the difference to the small system sizes 
which were accessible at that time (1997). No estimate 
for the critical point could be obtained from their work 4 . 
The data from Savenko and Dijkstra 2 - does not lie within 
the range of this graph. Presumably there is an error in 
reference [15] 30 . Therefore we cannot compare to this 
work in detail. 

The results presented here are of high accuracy: the er- 
rorbars are smaller than the symbols. The main sources 
of error are finite-size effects and insufficient sampling of 
the grand potential hypersurface perpendicular to the re- 
action coordinate (i. e., the packing fraction of rods). In 
order to check for the quality of sampling, we repeated 
our simulations for several different chemical potentials 
and different values of accuracy thresholds in the succes- 
sive umbrella sampling. From this we estimate the error 
on the coexistence volume fractions to be Arj r = ±0.002. 

Finite-size effects enter in two ways: 1) At coexistence 
the system forms slab configurations with two interfaces 
(due to the periodic boundaries). If L z is too small, the 
interfaces interact. Away from the critical point, we ruled 
out this contribution by increasing L z such that a plateau 
appeared between the peaks in P{rj r ). 2) P{rj r ) depends 
on the lateral box size, because the spectrum of capillary 
waves on the interfaces is cut off for wavelengths larger 
than the box. This effect is negligible far away from the 
critical point, too, because the interfacial tension there 
is large and the effects of capillary waves are very weak. 
However, close to the critical point a finite-size scaling 
analysis becomes necessary. In fact, the two effects are 
then just two aspects of the same property, the diverg- 
ing correlation length. We discuss this issue further in 
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FIG. 6: Cumulant ratio U/t as a function of the sphere fugacity 
close to the critical point for box lengths L x — 9D (crosses), 
12D (stars), 14D (squares) and 18D (triangles). Spherocylin- 
ders aspect ratio is L/D = 3. The inset shows a magnified 
plot of the region near the intersection point. 
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FIG. 7: Order parameter (r] l r — 77,?) (circles and diamonds) 
and interfacial tension 7 (triangles and squares) versus "rela- 
tive distance from the critical point" z s /zl — 1. Spherocylin- 
ders aspect ratios are L/D — 3 (triangles and circles) and 5 
(squares and diamonds). The solid lines correspond to the 
Ising power-law-behaviour with exponents Rising = 0.326 and 
2Hsing = 1-26. 



subsection II V A 21 

Close to the isotropic-nematic transition, the simula- 
tions cannot be equilibrated properly with the methods 
used here. For accurate grand-canonical simulations of 
the IN transition more advanced biasing techniques are 
necessary^. A study of this part of the phase diagram 
is beyond the scope of the present work, however. 



2. Finite-size scaling 

As explained above, P(rj r ) depends on system size. 
This is, however, not a drawback but an advantage, be- 
cause the finite-size effects can be used in order to locate 
the critical point of isotropic-isotropic demixing. One 
can construct quantities which are independent of sys- 
tem size at the critical point. One possible choic o 2 ' 13 of 
such a quantity is the cumulant ratio U4 — (m 4 )/(m 2 ) 2 
with m = T] r — (rj r ). When U4 is plotted versus z s close 
to the critical point for different system sizes, the inter- 
section point marks the critical point. Figure [S] shows 
that the critical sphere fugacity for spherocylinders with 
aspect ratio L/D = 3 is z c s = 1.109 ± 0.001. 

Figure[7]shows the difference in packing fraction r] l r —r]^ 
versus the "relative distance from the critical point" 
z s / z c s — 1 (circles for L/D = 3 and diamonds for L/D = 
5). The upper line in the graph is proportional to 
{z s /z c s — l)A=inB, where /?i s ing = 0.326 is the critical expo- 
nent of the order parameter in the Ising model. Clearly, 
the critical point is of the Ising universality class. 



3. Interfacial tension 

In the transition region where the gas phase transforms 
into the liquid phase and vice-versa, a grand potential 



barrier AVt needs to be crossed (indicated in figure [2]). 
Aft is related to the interfacial tension via 
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(24) 



where L 2 is the area of the interface and the factor 1/2 
accounts for the two interfaces, which are present due to 
the periodic boundary conditions 3 . Figure [7] shows values 
of the interfacial tension as a function of the "relative 
distance from the critical point" z s / ' z c s — 1 (triangles for 
L/D = 3 and squares for L/D = 5). The lower line in 
the graph indicates (z s /zf — l) 2l/Ising , where 2fi s j ng = 1.26 
is the critical exponent for the interfacial tension in the 
3d Ising model. 



B. Confinement 

1. Phase diagram 

Now we consider the behaviour of the mixture confined 
between two hard walls at distance d. Figures [5] and [S] 
show the phase diagrams for rod aspect ratios L/D = 3 
and 5. The distance between the walls is d/D = 3L/D. 
Demixing in confinement occurs at larger sphere fugac- 
ities than in the bulk. Also the chemical potential of 
the rods at coexistence is higher than in the bulk. The 
gas- like phase is shifted to larger rod volume fractions. 
The amount of spheres in the system is smaller than in 
the bulk at the same fugacity. The large increase of the 
concentration of rods in the gas-like phase distinguishes 
this system clearly from the behaviour of the Asakura- 
Oosawa-Vrij model in confinement. 

We also performed a finite size scaling analysis in con- 
finement to see how the critical point is shifted. Table U 
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FIG. 8: Phase diagram for a mixture of spherocylinders with 
aspect ratio L/D = 3 and spheres of diameter D between 
two hard walls at distance d/D = 3L/D (squares) compared 
to the bulk values (triangles). The filled symbols mark the 
critical points. The inset shows the phase diagram in the 
(z s , ?) r )-plane. The symbols show "raw data" for one box 
dimension L x = 12D only, and thus the curves marking the 
peaks of P(rj r ) do not join at the critical points ("finite size 
tails''-^). 
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FIG. 9: Phase diagram for a mixture of spherocylinders with 
aspect ratio L/D — 5 and spheres of diameter D between two 
hard wall at distance d/D = 3L/D (squares) compared to the 
bulk values (triangles). The filled symbols mark the critical 
points. The inset shows the phase diagram in the (z 3 ,T)r)- 
plane. The symbols show "raw data" for one box dimension 
L x = 20D only, and thus the curves marking the peaks of 
P(r) r ) do not join at the critical points ("finite size tails"—). 



L/D 


c 

*S,C.OBtf 


Vr,conf 


^s.bulk 


c 

Vr.bulk 


3 


1.26 ±0.01 


0.168 ±0.002 


1.109 ±0.001 


0.131 ± 0.002 


5 


1.05 ±0.01 


0.172 ±0.002 


0.981 ±0.001 


0.124 ± 0.002 



TABLE I: Critical point in the bulk and in confinement to 
slit-pore of width d/D — 3L/D. 




FIG. 10: Cumulant ratio U4 as a function of the sphere fu- 
gacity close to the critical point for box lengths L x = 12D 
(crosses), 15D (stars), 1873 (squares) and 21D (triangles) be- 
tween two hard walls at distance d/D = 3L/D. Aspect ratio 
L/D = 3. The inset shows a magnified plot of the region near 
the intersection point. 



The shift in rfc is caused by the wall-induced layer- 
ing of the rods. Because of orientational ordering, this 
effect is much stronger than e. g. the shift in the Asakura- 
Oosawa-Vrij model, which is due to positional ordering 25 . 
Close to the wall the rods are preferably oriented paral- 
lel to the wall ("parallel anchoring"). Therefore, their 
concentration is much higher than in the isotropic bulk. 

Obviously, it would be very interesting to investigate 
how the critical point depends on the wall separation d. 
This is unfortunately currently too demanding computa- 
tionally. One week of CPU time on a Pentium 4, 2.60GHz 
was needed to compute one density distribution in the 
(3L x 3L x 3L)-box, and one month of CPU time for the 
(5L x 5L x 5L)-box. Hence, a systematic study of the 
crossover from 3d to 2d in this system would require very 
large computational effort. 



lists the results in comparison to the bulk results. 

The relative shift in ?y£ , which we observe for rods with 
aspect ratio L/D = 5, is larger than the one of rods with 
aspect ratio L/D = 3, though the relative shift in z% 
behaves the other way around. The second effect is due 
to the ordering of the spheres close to the walls. Thus, 
its relative decrease on increase of the rods' aspect ratio 
is plausible, since, in terms of the sphere diameter, the 
distance between the walls increases. 



2. Order parameter profiles 

The walls change the structural properties of the gas- 
and liquid-like phase of rods. Here we show data from 
simulations in the A^UT-ensemble, where the volume of 
the system as well as the numbers of rods and spheres 
are fixed. The number of particles was chosen to match 
the coexistence values determined in the grand canonical 
ensemble. Although in principle finite size effects are dif- 
ferent in the canonical and grand canonical ensemble 13 , 
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FIG. 11: Density distribution of rods (triangles) and spheres 
(squares) normalised by their overall densities in the gas-like 
phase between two hard walls at distance d/D — 3L/D. 
Aspect ratio L/D = 5. Corresponding sphere fugacity is 

Z S = 1.1. 



FIG. 13: Nematic order S (triangles) and biaxiality £ 
(squares) parameter distributions in the gas-like phase be- 
tween two hard walls at distance d/D = 3L/D. Aspect ratio 
L/D = 5. Corresponding sphere fugacity is z s = 1.1. 
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FIG. 12: Density distribution of rods (triangles) and spheres 
(squares) normalised by their overall densities in the liquid- 
like phase between two hard walls at distance d/D = 3L/D. 
Aspect ratio L/D = 5. Corresponding sphere fugacity is z s — 
1.1. 
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FIG. 14: Nematic order S (triangles) and biaxiality £ 
(squares) parameter distributions in the liquid-like phase be- 
tween two hard walls at distance d/D = 3L/D. Aspect ratio 
L/D = 5. Corresponding sphere fugacity is z s = 1.1. 



far enough away from the critical point this difference 
can safely be neglected. 

To study the anchoring effects of the walls, we define 
the nematic order parameter S and the biaxiality pa- 
rameter £. S is the largest absolute eigenvalue^ of the 
matrix 

Q = y (3uV - I) , (25) 
' i— o 

where u l is a unit vector in the direction of the orien- 
tation of the rod i and I is the identity matrix. We 
divide the space between the walls into thin slices, so 
N r is the number of rods in such a slice at the distance 
z/D from the middle of the simulation box. (Note that 
many authors use the largest eigenvalue instead of the 



eigenvalue with the largest absolute value, which leads 
to different results in the case of uniaxial surface order- 
ing!) S indicates if there is a preferred direction in the 
system and how strongly the rods are oriented with re- 
spect to it. The eigenvector to this eigenvalue is called 
director. If S is zero, the phase is completely isotropic. If 
S is unity, all rods are aligned parallel to the director. If 
S is negative, they lie perpendicular to the director. The 
biaxiality measure £ is half of the difference of the other 
two eigenvalues of the matrix Q. It shows whether there 
is another preferred direction in the plane perpendicular 
to the director. 

Figures [TT] and [T2] show the density distributions of 
rods with aspect ratio L/D = 5 and spheres between 
walls at a distance d/D = 3L/D. Figure [Tl] shows 
the gas-like phase, figure [T2l the corresponding liquid-like 
phase at coexistence. The overall densities are approxi- 
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mately the positions of the peaks of the probability dis- 
tribution P(j] r ) from the grand canonical simulation at 
sphere fugacity z a = 1.1. 

Figures [T3] and [TJ] show the corresponding profiles of 
the nematic order parameter 5" and of the biaxiality pa- 
rameter £ in the gas- and liquid-like phases respectively. 

The positional as well as the orientational order of rods 
are clearly visible in the liquid- as well as in the gas-like 
phase. The range of the induced effects is of the order 
of the rod length, which is short in comparison to the 
chosen distance between the walls. In the middle of the 
system the order parameters reach their bulk values. The 
spheres are pushed away from the walls. This effect is 
also of the order of the rod length in both phases. 

V. DISCUSSION AND SUMMARY 

We have presented simulation results on the phase di- 
agram of mixtures of hard sphcrocylinders and "penetra- 
ble hard" spheres in the bulk and in confinement. We 
hope that these results are useful for experimental inves- 
tigations with suspensions of viruses and polymers. 

We have studied isotropic-isotropic demixing by sim- 
ulations in the grand canonical ensemble. In order to 
access states of high free energy we used the successive 
umbrella sampling method. The resulting phase bound- 
aries were compared to free volume theory. We extracted 
the critical point from an analysis of the cumulants of the 



1 T. Biben, P. Bladon, and D. Frenkel. Depletion effects 
in binary hard-sphere fluids. J. Phys.: Condens. Matter, 
8:10799, 1996. 

2 K. Binder. Finite size scaling analysis of ising model block 
distribution funktions. Z. Phys. B - Condensed Matter, 
43:119, 1981. 

3 K. Binder. Monte carlo calculation of the surface tension 
for two- and three-dimensional lattice-gas models. Phys. 
Rev. A, 25:1699, 1982. 

4 P. G. Bolhuis, A. Stroobants, D. Frenkel, and H. N. W. 
Lekkerkerker. Numerical study of the phase behaviour 
of rodlike colloids with attractive interactions. J. Chem. 
Phys., 107:1551, 1997. 

5 Y. L. Chen and K. S. Schweizer. Depletion interaction in 
suspensions of spheres and rod-polymers. J. Chem. Phys., 
117:1351, 2002. 

6 Y. L. Chen and K. S. Schweizer. Liquid-state theory of 
structure, thermodynamics, and phase separation in sus- 
pensions of rod polymers and hard spheres. J. Phys. Chem. 
B, 108:6687, 2004. 

7 M.A. Cotter. Hard spherocylinders in an anisotropic mean 
field: A simple model for a nematic liquid crystal. J. Chem. 
Phys., 66:1098, 1977. 

8 A. Cuetos, B. Martinez-Haya, S. Lago, and L.F. Rull. Use 
of parsons-lee and onsager theories to predict nematic and 
demixing behaviour in binary mixtures of hard rods and 
hard spheres. Phys. Rev. E, 75:061701, 2007. 



order parameter distribution. Free volume theory works 
well far away from the critical point, but, as expected, 
underestimates the concentrations at the critical point. 

In the bulk the system is very similar to the Asakura- 
Oosawa-Vrij-model. In particular, we showed that its 
behaviour on approach to the critical point falls into the 
Ising universality class. In confinement, however, the ori- 
entational degrees of freedom play a role. As the rods 
anchor parallel to the wall, the gas-like branch of the 
coexistence region moves to higher colloid (rod) volume 
fractions than in the Asakura-Oosawa-Vrij model. And 
the walls induce a much larger shift in the critical colloid 
(rod) volume fraction than they do for spherical colloids. 



Acknowledgments 

We would like to thank Jiirgen Horbach, Richard Vink 
and Peter Virnau for helpful suggestions. This work was 
part of the priority program SFB Tr6 (project D5) of 
the German Research Association (DFG). It was par- 
tially funded by the DFG Emmy-Noether-Program, the 
MWFZ Mainz. It has also been supported by the Eu- 
ropean Comission under the 6th Framework Program 
through integrating and strengthening the European Re- 
search Area. Contract: SoftComp VP-06/109. We 
thank the Forschungszentrum Julich for CPU time on 
the JUMP. 



9 Z. Dogic and S. Fraden. Ordered phases of filamentous 
viruses. Current Opinion in Colloid and Interface Sciences, 
11:47, 2006. 

10 T. T. Hebert. Phytopathology, 53:362, 1963. 

11 J. K. G.Dhont, M. P. Lettinga, Z. Dogic, T. A. J. Lenstra, 
H. Wang, S. Rathgeber, P. Carletto, L. Willner, H. Friel- 
inghaus, and P. Lindner. Shear-banding and microstruc- 
ture of colloids in shear flow. Faraday Discussions, 123:157, 
2003. 

12 K. G. Kang, M. P. Lettinga, Z. Dogic, and J. K. G.Dhont. 
Vorticity banding in rodlike virus suspensions. Phys. Rev. 
E, 74:026307, 2006. 

13 D. P. Landau and K. Binder. A Guide to Monte Carlo 
Simulations in Statistical Physics. Cambridge University 
Press, Cambridge, 2000. 

14 H. N. W. Lekkerkerker and A. Stroobants. Phase behaviour 
of rod-like colloid + flexible polymer mixtures. Nuovo Ci- 
mento D, 16:949, 1994. 

15 T. A. J. Lenstra, Z. Dogic, and J. K. G. Dhont. Shear- 
induced displacement of isotropic-nematic spinodals. J. 
Chem. Phys., 114:10151, 2001. 

16 M. P. Lettinga and J. K. G. Dhont. Non-equilibrium phase 
behaviour of rod-like viruses under shear flow. J. Phys. 
Cond. Matter, 16:S3929, 2004. 

17 W. Li and H. R. Ma. Depletion interactions between two 
spherocylinders. Eur. Phys. J. E., 16:225, 2005. 

18 R. J. Low. Measuring order and biaxiality. Eur. J. Phys., 



10 



23:111, 2002. 

F. Oosawa and S. Asakura. J. Chem. Phys., 22:1255, 1954. 
S. V. Savenko and M. Dijkstra. Phase behavior of a sus- 
pension of colloidal hard rods and nonadsorbing polymer. 
J. Chem. Phys., 124:234902, 2006. 

T. Schilling, S. Jungblut, and Mark A. Miller. Depletion- 
induced percolation in networks of nanorods. Phys. Rev. 
Lett, 98:108303, 2007. 

R. Tuinier, T. Taniguchi, and H. H. Wensink. Phase be- 
havior of a suspension of hard spherocylinders plus ideal 
polymer chains. Phys. Rev. E, page in press, 2007. 
R. van Roij. The isotropic and nematic liquid crystal phase 
of colloidal rods. Eur. J. Phys., 26:S57, 2005. 
R. Vink, S. Wolfsheimer, and T. Schilling. Isotropic- 
nematic interfacial tension of hard and soft rods: appli- 
cation of advanced grand canonical biased sampling tech- 
niques. J. Chem. Phys., 123:074901, 2005. 
R. L. C. Vink, K. Binder, and J. Horbach. Critical behavior 



of a colloid-polymer mixture confined between walls. Phys. 
Rev E, 73:056118, 2006. 

R. L. C. Vink and J. Horbach. The fluid-fluid interface 
in a model colloid-polymer mixture: Application of grand 
canonical monte carlo to asymmetric binary mixtures. J. 
Chem. Phys., 121:3253, 2004. 

R. L. C. Vink and T. Schilling. Interfacial tension of the 
isotropic-nematic interface in suspensions of soft sphero- 
cylinders. Phys. Rev. E, 71:051716, 2005. 
P. Virnau and M. Miiller. Calculation of free energy 
through successive umbrella sampling. J. Chem. Phys., 
120:10925, 2004. 

A. Vrij. Pure Appl. Chem., 48:471, 1976. doendnoteend- 
note30Dr. S. Savenko, private communication to Dr. R. 
Vink. 

Dr. S. Savenko, private communication to Dr. R. Vink. 



